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Abstract 

We review recent results from studies of the dynamics of classical Yang-Mills fields on a lattice. 
We discuss the numerical techniques employed in solving the classical lattice Yang-Mills equations 
in real time, and present results exhibiting the universal chaotic behavior of nonabelian gauge 
theories. The complete spectrum of Lyapunov exponents is determined for the gauge group SU(2). 
We survey results obtained for the SU(3) gauge theory and other nonlinear field theories. We also 
discuss the relevance of these results to the problem of thermalization in gauge theories. 
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I. INTRODUCTION 



Knowledge of the microscopic mechanisms responsible for the local equilibration of energy 
and momentum carried by nonabelian gauge fields is important for our understanding of 
non-equilibrium processes occurring in the very early universe and in relativistic nuclear 
collisions. Prime examples for such processes are baryogenesis during the electroweak phase 
transition, the creation of primordial fluctuations in the density of galaxies in cosmology, 
and the formation of a quark-gluon plasma in heavy-ion collisions. 

Whereas transport and equilibration processes have been extensively investigated in the 
framework of perturbative quantum field theory, rigorous non-perturbative studies of non- 
abelian gauge theories have been limited to systems at thermal equilibrium. We here review 
recent numerical studies of real-time evolution in the classical limit of lattice gauge theories. 
We will demonstrate that such an analysis can provide valuable insight into the dynamical 
properties of nonabelian gauge theories at high excitation energies. 

This paper is organized as follows. In section 2, we briefly discuss why a classical con- 
sideration of gauge fields is relevant, and we review the early studies of evidence for chaotic 
behavior of Yang-Mills fields. Section 3 is devoted to the results of our lattice study on 
SU(2) gauge theory. After some general considerations, we introduce the lattice formalism 
of gauge fields, and we describe the numerical techniques used to solve the field equations. 
Then we show the exponential divergence of two nearby trajectories, from which we can 
conclude that the system is chaotic, and extract the maximal Lyapunov exponent. In the 
last two subsections of section 3, we show how one can obtain the whole Lyapunov spectrum. 
After the experience with SU(2), results of several other theories are reviewed in section 4. 
These include compact U(l) and SU(3) gauge theory, massless Higgs fields, massive SU(2) 
vector fields, and finally the coupled SU(2) gauge-Higgs system. In section 5, we discuss 
one application of chaoticity in nonabelian gauge theory, namely, the thermalization pro- 
cess of highly excited gauge fields. We estimate the time needed to thermalize the system 
from the Kolmogorov-Sinai entropy of the gauge field. This time agrees with the results 
from thermal perturbation theory, supporting the view that the thermalization of the long 
wavelength modes is basically a classical process. Finally, in section 6, we point out some 
possible avenues of future work in this field, both in regard to the method itself and to 
physical applications. 
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II. GAUGE FIELDS IN THE CLASSICAL LIMIT 



A. General Considerations 

Let us start by discussing the classical limit of a simple system with a single degree of 
freedom which is described by a Hamiltonian 

H(p,x) = ~p 2 + V(x), (1) 

where V(x) is some well defined potential. The quantum evolution is defined either by 
operator equations for p and x in Heisenberg picture, 

x = p, 

P = d x V(x), (2) 
or by a Schrodinger equation for the wave function, 

ihdt$ = H®. (3) 

Both are appropriate to describe the behavior of a pure state. To admit mixed states it is 
necessary to use the quantum Liouville equation 

ihd t p=[H,p] (4) 

for the density matrix p. For a system in a pure state, p = |$)($|, the Liouville equation is 
reduced to the Schrodinger equation. 

Quantum mechanics and classical mechanics are related by the correspondence principle, 
which states that at large quantum numbers or small % quantum mechanics has classical 
mechanics as a limit. But it is not trivial to actually construct the correspondence, especially 
when the system considered is chaotic in the classical limit. Quantum chaos is still far from 
being a well understood concept 0. 

In the Heisenberg picture the classical limit is attained by treating the operator equations 
as equations of real numbers, i.e. by neglecting the non-commutivity of x and p which is 
of the order of h. The classical limit is equivalent to the limit h — > 0. The advantage of 
working in this picture is that the classical limit is directly reached by the above simple 
prescription. The disadvantage is that it is difficult to estimate the accuracy of the classical 
approach. 

In the Schrodinger picture, it is useful to consider wave functions which have minimum 
spread in both momentum and space representation. When the system is highly excited, i.e. 
the accessible phase space is much larger than the volume of the wave packet, it describes 
a state for which both x and p have "sharp" values. The classical equations of x and p are 
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obtained by calculating time derivatives of expectation values of x and p under such a wave 
function evolving according to the Schrodinger equation and neglecting the corrections due 
to the finite width of the wave packet. These corrections can be shown to be of the order of Ti 
and hence can be neglected for a highly excited system. The problem in this representation 
is that as the system evolves with time, the minimum uncertainty may not be maintained 
due to the influence of interactions. The width of the wave packet may start to increase and 
the classical equations may become invalid. This problem is especially serious for classically 
chaotic systems. An exact quantum calculation has been done for the "Arnold Cat" system 
which is chaotic in the classical limit. It was shown that an initially sharp Gaussian wave 
packet quickly dissolves into a diffuse state which is anything but a Gaussian jg. Hence 
we must be especially cautious when we talk about a localized wave packet in a chaotic 
quantum system. 

Another way to compare classical and quantum mechanics is to work with the Wigner 
function W(p,x) which is defined as 

/+oo 1 1 

^ dye- ip y/ h (x + -y\p\x--y). (5) 

From the Liouville equation we can derive the time evolution of W(p,x), 
dW 

— = (£ c + £ q )W (6) 

where C c is the classical Liouville operator which describes the classical dynamics in phase 
space language, 

C c = (d p H)d x - (d x H)d p , (7) 
while C q comprises the quantum corrections 

c q = f A {dlv)&l-^{dlv)dl + ... (8) 

which are of successively higher order in %. In principle, in this representation the quantum 
corrections to the classical calculation can be computed order by order in h. The rapid 
dissolution of localized wave packets for chaotic systems has also been observed in the Wigner 
representation 0], and some attempts have been made to include quantum corrections by 
means of a stochastic process |5|. The Wigner function approach has been invoked to 
derive transport equations for gauge theories in the mean- field approximation 6], but to our 
knowledge it has not yet been applied to gauge theories on a lattice. 
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B. Classical Limit of a Gauge Theory 



Now let us turn our attention to nonabelian gauge fields, where we are dealing with a 
system of infinitely many degrees of freedom. What does it mean in this case that the 
system is highly excited? Certainly not every degree of freedom can be highly excited, 
because this would require an infinite amount of energy. To this end it is better to look at 
the system in Fock space. The system is excited by creating particles with different energy 
E and momentum hk. Now suppose the gauge field is excited to a temperature T. The 
Bose distribution function implies that the long wavelength particles are more copiously 
excited. For E <C T the Bose distribution (e E ^ T — merges into classical distribution 
function T/E, whereas the "hard" particles of short wavelength are rarely excited. So, 
generally, the long wavelength modes can be treated classically and the short wavelength 
modes retain their quantum statistics. With increasing T, or decreasing H, more and more 
modes approach the classical limit. If, at certain temperature, a physical quantity is only 
related to long wavelength modes, then we expect the classical calculation to be adequate 
and the quantum corrections to be small. 

It is here where the lattice regularization of the gauge field plays an essential role. While 
being originally advocated Q as gauge-invariant cut-off of the ultraviolet divergencies of 
the quantum field theory, it can assume a new role in the definition of the classical high- 
temperature limit of the quantum field theory. In the lattice formulation all modes with 
wavelength shorter than the lattice spacing a are eliminated, and we are left only with the 
infrared modes. Hence the lattice regulated gauge theory generally goes to the classical limit 
at high temperature T 3> h/a or vanishing %. The validity of the classical calculation may 
depend on the nature of the quantity we are interested in. For example we will show later 
that the damping rate for a gluon at rest, as well as the rate of thermalization, are essentially 
classical quantities, because they are induced by the interaction of long wavelength modes. 
On the other hand, the screening length of static electric gauge fields is controlled by short 
wavelength modes , and hence is at best a semi-classical quantity. 

The second problem is how to derive the classical field equation from the quantum field 
equation. Does the classical equation make any sense? This question is related to the 
confinement problem. We can see this as follows. In the Heisenberg picture we can start 
with the operator equation 

= 0, (9) 

where D is the covariant derivative and F is the field operator. The classical equation is 
obtained by treating F and A as numbers rather than operators. We call F(x) and A(x) 
classical field configurations. Suppose now we want to work in the Schrodinger picture. In 
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an unbroken nonabelian gauge theory, for example in QCD, we observe color confinement 
so that only color singlet states are physically realizable. If we calculate the expectation 
value of a (color-octet) electric field component E a (x) in a singlet state, the result vanishes 
according to the Wigner-Eckart theorem. This means that a classical configuration does 
not correspond trivially to a physical state with minimum uncertainty. The reason is that a 
classical gauge configuration is only gauge covariant while a singlet state is gauge invariant. 
But this does not mean that there is no relation between them. We can start from a classical 
configuration and gauge rotate it to obtain other configurations. We then superimpose these 
to form a singlet state which corresponds to a physcial state. This procedure corresponds to 
the Peierls-Yoccoz projection method [8] for wavefunctions with good symmetry properties. 
We can reverse the process and decompose a highly excited color singlet state into some 
"nearly classical" configurations. By studying these configurations classically, we can gain 
insight into the evolution of the original physical state. 



C. Reasons for Investigating the Classical Limit 

The most commonly used approximation method in quantum field theory is perturbation 
theory. But we know that in nonabelian gauge theory there exist some fundamental non- 
perturbative effects, such as color confinement and topological quantum numbers. These 
effects are beyond the reach of perturbation theory; non-perturbative methods are required 
to understand them. The investigation of gauge fields in the classical limit provides one 
such approach. Studies of the classical field equations have led to some very interesting 
non-perturbative results, such as monopole solutions and instantons js, 10]. These classical 
solutions prove to be vital to the understanding of the corresponding quantum physics. On 
the other hand, these solutions are not general integrals of the classical non-linear gauge 
field equations, but exploit special symmetries of the nonabelian gauge theory. In any case, 
the known classical solutions of the Yang-Mills equations are by no means exhaustive, and 



the equations have been shown to be nonintegrable in general ll|. Thus we do not have 
a complete understanding even for the classical field equations. It is our hope to learn 
something about the general behavior of nonabelian gauge fields by numerically integrating 
the classical equations of motion. 

Secondly, under some extreme conditions, e.g. at high temperature, the quantum field 
reaches its classical limit at least for some observables. Quantities that are calculable in the 
classical limit can be identified in the high-temperature expansion of perturbation theory as 
those that exhibit a leading term proportional to g 2 T. Namely, if g is the coupling constant 
of the classical gauge theory, related to the standard dimensionless gauge coupling constant 
a = g 2 h/4:7T, the term g 2 T has the dimension of an inverse length or time and survives in 
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the limit % — > 0. Quantities with this leading behavior are, e.g., the gluon damping rate [12 1 
and perhaps the inverse screening length of static magnetic gauge fields 13J, [l^ . If we are 
interested in these quantities, the classical approach can provide us with a practical method 
of calculating these physical observables. 



D. Chaotic Dynamics of Yang-Mills Fields 

The first evidence that Yang-Mills fields exhibit chaotic dynamics was found a decade 
ago by Matinyan, Savvidy, and others [ill, Q, 13, > when they studied the dynamics of 
spatially constant potentials in the SU(2) gauge theory. In contrast to electrodynamics, such 
potentials are not always gauge equivalent to the trivial vacuum, due to the non- commutative 
nature of gauge transformations in nonabelian gauge theories. The original motivation for 
these studies was the desire to show that the Yang-Mills equations form a non-integrable 
dynamical system; it obviously suffices to prove this assertion in a limiting case. However 
one can also regard constant potentials as the relevant degrees of freedom surviving in the 
infrared limit. Indeed, Liischer has shown that the Hamiltonian (1) appears as lowest order 
term in the effective action for the Yang-Mills theory on a three-dimensional torus, i.e. in a 
cubic box with periodic boundary conditions [3]. 

Choosing the temporal gauge A% — (a = 1, 2, 3), and assuming that the vector potentials 
A a (t) are functions of time only, the dynamics is governed by the Hamiltonian 

#™ = E l^y + is 2 E( Aa >< A ') 2 , (io) 

a a, ft 

where g is the gauge coupling constant. It is not hard to show that the system allows for 
seven integrals of motion, corresponding to energy, angular momentum, and color charge 
conservation. In fact, the Hamiltonian can be reduced to the form 

H YM ~ \(x 2 + y 2 + z 2 ) + \g 2 {x 2 y 2 + y 2 z 2 + z 2 x 2 ) + ... (11) 

where the dots indicate terms describing quasi-rotational degrees of freedom. The nontrivial 
dynamical aspects are all contained in the three variables x(t), y(t), z(t). Note that the 
coupling constant g can be eliminated by rescaling the time coordinate. It is therefore useful 
to introduce an additional term into the Hamiltonian which breaks this scale invariance, e.g., 
a harmonic potential: 

H YMH = H Ym + \g 2 v 2 {x 2 + y 2 + z 2 ). (12) 
The dynamical properties of this Hamiltonian are controlled by the dimensionless parameter 
r = \{gvY/g 2 E 1 (13) 
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where E is the energy density. Numerical studies of Poincare surfaces of section of trajec- 
tories revealed ll| that the motion governed by the Hamiltonian ()12j) is regular for values 
r ^> 1, becoming partially chaotic as r falls below 1, and strongly chaotic for r — > 0, as 
shown in Figure 1. The mechanism leading to divergence of nearby trajectories is similar 
to a classical billiard: the equipotential boundary of classical motion for a given energy E 
has negative curvature. A singular-point analysis of the system (|T2*|) has been performed by 



Steeb et al. [20]. It is not known whether the Hamiltonian describes a true K-system, 
this conjecture has been refuted for the two-dimensional analogue |2l| 

Space-dependent solutions of the Yang-Mills equation have been studied in the case of 
spherical symmetry for solutions of the form J22, 23, 24] 

which include the so-called Wu-Yang monopole (( 
the nonlinear wave equation 



(14) 

0). The radial function <f)(r,t) satisfies 



(15) 



which has been shown to exhibit the rapid energy sharing between Fourier components 
which is characteristic of chaotic systems 2^0- Equation (fTo]) has also been shown to be 
non-integrable by the method of Painleve analysis 24 1 . 

Recently, solutions of the Yang-Mills equations in two spatial dimensions have been stud- 
ied numerically, subject to the assumption that the potentials depend only on one spatial 
coordinate and on time 25j . Again mode-sharing of the energy was observed, and there are 
indications that the spatial potential functions evolve into a fractal pattern. 



III. CHAOS IN SU(2) LATTICE GAUGE THEORY 
A. General Considerations 

Interesting as these results are, they leave two important questions unanswered: What 
are the dynamical properties of the full (3+l)-dimensional classical Yang-Mills field? What 
is the physical significance of chaotic dynamics of the classical field theory? The goal of our 
investigation was to provide at least partial answers to these questions. Our approach [2^ 
deviates in two important aspects from earlier studies: 

(a) Since the spatial coordinates have to be discretized for numerical purposes, it is 
convenient to formulate the SU(2)-gauge theory in terms of matrix- valued link variables on 
a N 3 cubic lattice 0,Ezj] : 

U Xti = exp (-ligaA^x)^) . (16) 
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Here r a are the Pauli matrices, a is the elementary lattice spacing, and (x, i) denotes the 
link from the lattice site x to the nearest neighbor in direction i, x + i. The link variables 
U x> i are explicitly gauge covariant, as opposed to the Yang-Mills potentials Af(x). Since the 
U x>i take values on the gauge group SU(2) rather than the group algebra, the magnetic field 
strength is bounded for a given lattice spacing. 

(b) Instead of studying the gauge field dynamics in the vicinity of arbitrarily selected 
configurations, we have investigated the dynamical behavior of random field configurations, 
corresponding to gauge fields selected from a microcanonical or canonical ensemble. This 
has the advantage that our field configurations are controlled by a single parameter, the 
temperature T or the average energy density e, which can be varied systematically. This 
approach also allows for the identification of some quantities calculated for the classical 
Yang-Mills theory with those obtained in the high-temperature limit of the quantum field 
theory. 



B. Lattice SU(2) Theory in Hamiltonian Formalism 



Our study is based on the Hamiltonian formulation of lattice SU(2)-gauge theory 28, 
governed by the Hamiltonian 

H = %I>(Ptfi«d + 4-Et 1 " |trC/^], (17) 

y x,i y a x,ij 

where a dot denotes the time derivative. Here the electric and magnetic fields have been 
expressed in terms of the SU(2) link variables U X) i(t) and the so-called plaquette operator 
U Xy ij which is the product of all four link variables on an elementary plaquette with corners 
(x, x + i, x + i + j, x + j): 

with U x -i = Ul_ Li . The links are directed and hence the plaquettes are oriented. In the 
continuum limit the plaquette variable U x> ij is related to the local magnetic field k 

U x ,ij = exp(-i^ga 2 e ijk B^ k T a ), (19) 
while the electric field on the lattice is given by 

Ki = -^^ atJ ,M,d- (20) 

In the classical limit, this Hamiltonian is scale invariant. To see this explicitly, we scale the 
time variable, s = t/a, obtaining 

y 2 Ha = £tr 6 -^) +431- ±tr U^}, (21) 

x,i \ / x,ij 
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where the right-hand side is parameter free. So the only parameter in the system is total 
energy or temperature. 

The classical equations of motion are derived from this Hamiltonian, making use of the 
lattice representation of the electric field components (|20j). It is useful to write the SU(2) 
matrices in the form 

rr ■ ( U - iU 3 , U 2 - IU X \ 

U = u - ir a u a = . . (22) 

I — u 2 — lUx, U + lUz J 

where the Uj are four real numbers, which can be thought of as components of a quaternion. 
The unit determinant implies that 

det U = ul + u\ + u\ + u\ = 1. (23) 

One easily verifies that the components of the quaternion satisfy the following differential 
equations of motion: 

(u U = Y a E >U ( 24 ) 

and 

2 



C \ -!L 
[U a ) x ,i ~ ^ 



IZ,(»»lr, ■ <""' ■ (25) 
These conserve the quaternion length, \\U\ \ = detU, because they satisfy 

UqUq + il a U a = 0. (26) 

Although it would be sufficient to update only the three fields u a during the time evolution 
because of the unit length constraint, it is computationally more efficient to also update the 
fourth component Mo, as well using eq. (J24|) rather than calculating uq by taking the square 
root of (1 — u a u a ). The three electric fields on each link are updated according to 

Ki = —2^ [^ a (U X!ij - Ul 3 )] (27) 
ag ■ z 

where the sum runs over all four plaquettes that are attached to the link (x,i). Here 
tr (^r a Q), for any quaternion Q, just reads off the component q a , and therefore does not 
require additional computational effort. 

We note here that the time evolution for the electric field conserves Gauss' law 

BfE% = 0, (28) 

which is an expression of charge conservation. Starting with any configuration (U x ,i, E x ^) 
that satisfies Gauss' law, the time evolution will not violate it. Precisely because of the 
chaotic nature of Yang-Mills dynamics it is, in general, impossible to integrate this set of 
equations analytically. We have to rely on numerical methods, which will be discussed in 
the following subsection. 
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C. Integrating the Equations of Motion 

1. Lattice Geometry and Boundary Conditions 

In all our simulations we discretize the gauge fields on a simple cubic iV 3 lattice with 
periodic boundary conditions. To calculate the oriented plaquette products more efficiently 
we use a linked-list approach. All the sites are numbered consecutively. The link numbers are 
in x, y and z order, and are arranged according to the site of origin. The plaquette numbers 
are chosen so that they are orthogonal to the links in x, y and z directions, respectively. 
We also build two additional arrays of numbers at the start of the simulation. One contains 
the link number of all the links that form a given plaquette, ordered in the way they appear 
in the directed plaquette product. The other list contains the number of plaquettes that 
contain a given link, ordered according to the position of the link in the plaquette. The 
periodic boundary conditions are reflected automatically in the above two link lists. 

2. SU(2) Representation 

There are several ways to represents a SU(2) matrix: by a complex 2x2 matrix, by a 
real quaternion, or in polar coordinates as a point on the 4-dimensional unit sphere which 
can be specified by 3 angles. The last representation has the lowest storage requirements, 
however, it involves time-consuming trigonometric conversions. 

In our simulations we chose the quaternion representation because it takes the least 
number of floating point operations to multiply two group elements. In the complex matrix 
representation one needs 32 multiplications and 24 additions compared to 16 floating point 
multiplications and 12 additions when taking the product of two quaternions. Moreover, 
the equations of motion can be written in a simpler form in quaternion representation. 

3. Numerical Integration 

The numerical task consists in propagating the gauge fields in time, by integrating the 
equations of motion. This can be achieved by a variety of numerical methods. Here we use 
the Runge-Kutta method with fourth-order accuracy, which is easy to implement, allows for 
adjustable time-step control, and is quite stable. 

The most CPU intense part of the simulation is the computation of the oriented product 
of the fields over the complement lattice. This complement lattice is defined by all oriented 
links contained in the elementary plaquette attached to a given link. The overall performance 
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of the code depends on several issues like field representation, fast access to the data on the 
neighboring links and plaquettes, etc. 

4- Code Verification and Accuracy 

Several integrals of motion can be used to verify the simulation code and to test the 
accuracy of the time integration. Since we are studying the Hamiltonian evolution of the 
gauge fields, the total energy of the system remains constant. With a typical value for the 
time step of the order of 0.01 dimensionless units, the total energy is conserved to better 
than 8 significant digits. 

Another conserved quantity is the length of each quaternion link variable for SU(2). For 
an interval less than one dimensionless time unit the conservation is better than 12 significant 
figures. However, due mainly to accumulation of cut-off and rounding errors the precision 
is deteriorating progressively. Therefore we choose to rescale the gauge fields after every 
time integration to maintain fixed length of each dynamical variable, since the subsequent 
integrations are very sensitive to the quaternion length preservation. This method permits 
to integrate the equations of motion with a larger time step. Later we will also study SU(3) 
gauge theory, where the corresponding criterion is that the determinant of the unitary matrix 
on each link shall be unit. In this respect, the SU(3) evolution is more stable. The deviation 
of the determinant from one is of order 10~ 8 even if the system evolves over a total time 
interval of T = 20 or 30. 

The validity of Gauss' law, eq. (|28|1 . is also an indication for accurate integration, which 
is extremely sensitive to all kinds of program errors and thus provides a valuable probe for 
code verification. In our calculations the color charge was always conserved to better than 
five significant digits. 

5. Performance 

The most time consuming part of the code is the calculation of oriented plaquette prod- 
ucts. However, this part of the code is fully vectorizable and runs at about 160 Mflops 
for SU(2) and 100 Mflops for SU(3) on a single Cray-Y-MP processor. Typically a single 
time-step integration of the set of equations of motion for N = 10 takes about 30 ms CPU 
time for SU(2) and 220 ms for SU(3) on the same processor. 
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D. Divergence of Trajectories 



In order to look for chaotic motion we will consider the evolution of infinitesimally sep- 
arated gauge field configurations. If we find exponentially diverging trajectories, we can 
identify the positive Lyapunov exponents and obtain the value of the Kolmogorov-Sinai en- 
tropy, i.e. the entropy growth rate, of the gauge field. To observe the exponential divergence 
of two trajectories U x ^{t) and U' x i {t) in the space of gauge field configurations we introduce 
the following gauge-invariant distance: 

DP^U'J = ^£|tr[4 !ii -tr[/; i .|, (29) 

ZiV P x,ij 

where N p = 3iV 3 is the total number of elementary plaquettes. In the continuum limit 
DP^U'J a ^ D{A'lXt}K ^ j d?x\B\xy - B{x)\ (30) 

i.e. D measures the average absolute local difference in the magnetic energy of two different 
gauge fields. We note a peculiar property of this distance measure, which is a natural 
consequence of the topology of the compact 9iV 3 -dimensional space of magnetic gauge field 
configurations on the lattice: For N ^> 1 almost all pairs of configurations have the same 
distance D. This is illustrated in Figure 2, where the distribution of distances of randomly 
chosen configurations from a fixed field configuration is shown for lattices of different size. 
For large lattices the distribution approaches a narrow Gaussian with a width of order N 3 / 2 . 
This property does not limit the usefulness of the metric (}2*9"j) as measure of the divergence 
of infinitesimally separate field configurations, but it causes the saturation of D at large 
times observed in the calculations (see figures below). 

Figure 3 shows the evolution of D(t) for initially neighboring gauge field configurations 
on a 20 3 lattice. We choose the reference configuration by randomly selecting link variables 
in such a way that the average energy per plaquette takes on the desired value This 
procedure is controlled by a parameter 5 which varies between and 1. The energy per 
plaquette grows like 5 2 for small S and saturates in the limit 8 — > 1. We then construct 
a neighboring configuration by perturbing each link element infinitesimally (see [26] for 
details). For values of 5 of order unity the distance D{t) starts to grow exponentially as 
D{t) = D exp(ht) almost immediately (see Figure 3a). The growth rate h decreases with S, 
and for 5^1 one observes an extended period during which the distance D(t) between two 
adjacent field configurations performs more or less regular oscillations before exponential 
growth finally sets in (see Figures 3b, 4). For very small values of 8, corresponding to very 
low energy density of the gauge field configurations, the exponential growth pattern of D(t) 
is modulated by low-frequency oscillations, which we attribute to the growing influence of 
non-leading Lyapunov exponents. 
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The initial latency period before the onset of exponential growth of D can be estimated as 
follows: We can write U' xi {t) = U x ^{t) + 5U x ^{t) where 5U X ^ approximately satisfies a system 
of linear differential equations with 18iV 3 eigenmodes and eigenfrequencies. On average every 
eigenmode will be excited with equal probability by our random choice of 8U Xt i(0). In order 
for the maximally unstable mode to outgrow the combined weight of all other modes we 
therefore have to wait a certain time to, statistically given by exp(hto) « (lSiV 3 )^, i.e. 



This agrees roughly with the observations, in particular, it explains why the onset of expo- 
nential growth is delayed for small values of the slope parameter h (see Figure 4). 

E. Maximal Lyapunov Exponents 

It is natural to identify the growth rate h with the maximal Lyapunov exponent Ao of 
the lattice gauge theory. This is confirmed by a careful analysis of the Lyapunov spectrum 
(see sections 3.6 and 3.7 below), and we will assume the equality h = Ao in the following. 
Extensive studies have shown that Ao is a universal function of the average energy per 
plaquette E, as shown in Figure 5. For values 5 > 0.15, the numerical determination of 
Xo(E) from D(t) is quite reliable, and the statistical and systematic errors are small. Figure 
5 demonstrates that Xq(E) is growing approximately linearly with energy. Using the property 
of scale invariance of the Hamiltonian, discussed previously in eq. (|21jl. one finds that the 
dimensionless product Aoa can only be a function of the combination g 2 Ea. Our numerical 
results show that this function is approximately linear: 



This scaling property has been verified numerically over a wide range of values for g and a 
(see Figure 5). We note that, according to (}32|) . Ao is independent of the lattice spacing a 
in the classical limit, where g does not depend on a. 

We have also studied the dependence of X (E) on the size of the lattice at fixed lattice 
spacing a. Figure 7 shows the evolution of the distance between adjacent field configurations 
for lattices of size ranging from N = 6 to iV = 28. The rapid convergence is obvious, the 
curves for N = 28 hardly deviate from those for N = 6, except for a decrease of fluctuations 
which is not visible in the figure because the curves almost coincide. We have not observed 
within statistical errors a systematic dependence of Xq(E) on N, for N > 6 and S > 0.15. 
The results obtained for N = 6 are shown in Figure 5 as solid squares. For smaller values 
of S we found that the exponential growth rate gradually decreases with growing N, which 
may explain why the lowest two points in Figure 5, obtained for N = 20, still lie above the 
straight line. 



t ~ ln(18iV 3 )/2/i. 



(31) 



A a « \g 2 Ea. 



(32) 
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F. Rescaling Method for Lyapunov Analysis 



The previous method for obtaining the largest Lyapunov exponent of a Hamiltonian 
system is straightforward, but it has two drawbacks. First, the exponential divergence of 
trajectories shows fluctuations, resulting in an uncertainty in the determination of the ex- 
ponential divergence rate. The second drawback of the method is that only the largest 
Lyapunov exponent can be obtained in this way, but not the complete Lyapunov spec- 
trum. We now discuss a method which can be used to determine Lyapunov exponents more 
precisely and yields the whole spectrum of Lyapunov exponents. 

This technique, which we call rescaling method here, is widely used in studying chaotic 
dynamical systems Suppose we want to calculate the two largest Lyapunov exponents 
of the system. We can randomly choose three initial points in phase space, to which we 
refer as zq(0), zi(0), £2(0) for convenience, with the condition that they are close to each 
other according to some appropriate distance measure. If the system is chaotic, or if the 
initial points are chosen inside the chaotic part of the phase space, the distances between 
the three trajectories Zi(t) evolved from the three initial points will diverge exponentially. 
Let us denote the separation vectors between the trajectories by 

di(t) = Zi (t) - z (t), (i = l,2) (33) 

and the absolute distance between Z{ and zq by Di = Since the available phase space 
volume is limited by the total available energy, Di(t) will saturate after a certain time. 
To avoid saturation, the following rescaling method is used. The fixed reference distance 
D = Di(0) is chosen in the beginning. The whole procedure consists of two steps. In the 
first step, the trajectories Zi(t) evolve according to the equations of motion for a period of 
time t . Then in the second step we rescale the separation vectors d\ and d>2 as follows. We 
hold the point z (t ) fixed, but scale the distance D 1 back to D Q by setting: 

z[(t ) = z (t ) + -^L-d^t,,). (34) 

The scaling factor is denoted by 

4 = D 1 (t )/D , (35) 

where k refers to the kth rescaling. For ^2(^0); we first orthogonalize 0^2(^0) against di(t ), 
then scale the orthogonalized vector d! 2 to the reference length Dq. The new scaling factor 
is denoted by s|. This procedure is iterated n times until the Lyapunov exponents, 

n i k 

^JS&E^. (36) 

k=l L o 
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converge. Here Ai is the largest Lyapunov exponent and A2 is the second largest one. By 
adding more and more trajectories, in principle, this method can be used to obtain the whole 
Lyapunov spectrum with arbitrary accuracy. The time needed to obtain a given Lyapunov 
exponent depends on how fast the procedure converges. 

At first glance, the applicability of this method to lattice gauge theory is questionable. 
The reason is that while the meaning of rescaling and orthogonalization is quite obvious in 
a Euclidean phase space, it is less clear in the case of a gauge field, where the phase space is 
curved as well as constrained by Gauss' law. The problem of curvature is relatively simple. 
One possible approach is to transform the link variable U Xji back to the vector potential A*, 
and to work in the phase space formed by electric fields and vector potentials, for which the 
geometry is Euclidean. In the case of SU(2), we have yet a simpler method. Each SU(2) 
group element is represented by a normalized quaternion. When performing orthogonaliza- 
tion and rescaling, we can simply treat all components of the quaternion as independent 
cartesian coordinates, because locally the metric of the curved space is Euclidean. 

Surprisingly, Gauss' law also does not pose a serious threat, for the following reason. 
When we rescale the separation vectors d i: we indeed violate Gauss' law. But if the distance 
Di is small, then the violation of Gauss' law is of second order in A- If we limit ourselves to 
sufficiently small Di(t ), then the violation of Gauss' law in each rescaling step is negligible. 
We next observe that the evolution of the system respects Gauss' law, so the violation 
does not increase with time. On the other hand, the next rescaling decreases the previous 
violation of Gauss' law by a (large) scale factor, so the violations do not accumulate. The 
same argument applies to the small changes in the choice of gauge induced by the linear 
rescaling procedure. 

We have used this method to study the SU(2) theory and have measured the largest two 
Lyapunov exponents. We indeed find that the violation of Gauss' law remains of the order 
of 1CT 6 . The result of a typical run for the Lyapunov exponents is shown in Figure 8 for 
a configuration with scaled energy g 2 Ea = 4.06. The solid line corresponds to the largest 
exponent and the dotted line to the second largest one. They converge at f « 100. Note 
the time scale of saturation of the distance D(t) in the case without rescaling is about 30 
at the same energy. The result obtained with our new, improved method, A = 0.667, is 
very close to, but slightly lower than, our previous result 4.06/6 = 0.677, where we did not 
use the rescaling technique. We note that the results for the Lyapunov exponents generally 
converge from above, i.e. the Lyapunov exponents are overestimated when the trajectories 
are not followed for sufficiently long time. We also observe that Ai is almost identical to A . 
The reason is that, as we are going to show next, there exists a whole Lyapunov spectrum 
which forms a continuous curve in the thermodynamic limit. 
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G. Lyapunov Spectrum 



The above method, in principle as well as in practice, can be extended to obtain the 
whole Lyapunov spectrum for small lattices. For each additional Lyapunov exponent, we 
need to integrate one more trajectory so as to form one additional linearly independent vector 
cfj. If we want to calculate vl Lyapunov exponents, we need to run ut, + 1 configurations 
simultaneously, from which we can form ul separation vectors and calculate the Ul largest 
Lyapunov exponents. Practically, the computational resources put limits on the possible size 
of the lattice for which we can obtain the whole Lyapunov spectrum. We have performed 
studies with l 3 , 2 3 and 3 3 lattices. Fortunately, we found that the spectrum starts to scale 
as early as size 3 3 , which permits us to extrapolate the result to the thermodynamic limit. 

We have calculated the complete spectrum for l 3 and 2 3 lattices. The results for the latter 
are shown in Figure 9. On a N 3 lattice the dimension of phase space is 9N 3 x 2, because 
there are three vector and three color directions at each site for magnetic and electric fields. 
This amounts to a total number of 144 Lyapunov exponents for N = 2. We can see that 
the spectrum is divided into three equal parts. The first one third of Lyapunov exponents 
are positive, while the second one third are all zero, and the last one third of exponents are 
the negative of the first one third. The vanishing Lyapunov exponents account just for all 
the degrees of freedom associated with static gauge transformations and with conservation 
of Gauss' law, each equal in number to three times the number of lattice sites. Thus our 
results confirm the general properties of Lyapunov spectra|34]. The results for a l 3 lattice 
are basically the same, but the Lyapunov spectrum consists only of 18 numbers. The value 
of the largest Lyapunov exponent is consistent with the result obtained earlier for the model 
of spatially constant Yang-Mills potentials (see section 2.4). 

In Figure 10 we show the scaling of the Lyapunov spectrum, where we compare the results 
from a 2 3 lattice and a 3 3 lattice. To save computation time, we have calculated only the 
positive Lyapunov exponents for the 3 3 lattice. In both cases, the initial configurations are 
chosen similarly. For presentation the Lyapunov exponents were scaled with respect to the 
largest exponent Ao in each case. The indices labeling the Lyapunov exponents are scaled 
to the total number of Lyapunov exponents, i.e. 144 for a 2 3 lattice and 486 for a 3 3 lattice. 
The two lines coincide nicely, exhibiting an early scaling behavior. 

The scaled Kolmogorov-Sinai entropy, i.e. the sum over all positive Lyapunov exponents, 

is: 

§^ 2 ' < 37 > 

where N 3 is the size of the lattice. We point out that at this small lattice size, the spectrum 
does not yet scale with energy. But we know from section 3.4 that this scaling appears at 
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a lattice of size 6 3 , where we found X /g 2 E ks |. If we conbine both scaling laws, we can 
estimate the Lyapunov spectrum in the thermodynamic limit, from which we obtain the 
Kolmogorov-Sinai entropy density 

*KS = ^EA,«i/ £ , (38) 

where e = 3E/a 3 denotes the energy density, and E is the average energy per elementary 
plaquette. Makinmg use of the thermodynamic relation oT = e + P = |e, we find that the 
characteristic entropy growth rate for SU(2) is given by 

^ " \^ = T2 9 ' 2r ' < 39 > 

er y a 12 

The above method is quite successful for SU(2) gauge theory, but it is not obvious how 
to apply it to SU(3) gauge fields. The reason is that this method relies on the quaternion 
representation, which is quite special for the group SU(2). We have also developed a more 
general method to obtain the Lyapunov spectrum, which can be used to study SU(3) gauge 
theory. Again the basic idea is simple. We want to directly consider the motion of vectors 
in the tangent space built upon the phase space {E xi ,U Xt i}. The tangent space of E Xji is 
simple, the vectors just being given by SE x i . We must be more careful specifying a vector 
in the tangent space to U X j. Here, to be consistent with our definition of the conjugate 
momenta E x>i as the left group generators, we define a vector b x ^ in the tangent space of 
U X j as 8U X j = ib X) iU Xy i. A vector in the complete tangent space is the direct product of b x ^ 
and 5E Xj i. The linear evolution equations for b x ^ and 8E X ^ can be derived and integrated 
along with the equations of motion in phase space. Initially choosing vl vectors, we can 
again obtain Lyapunov exponents. The results for SU(2) by this method is shown in 
Figure 11, in comparison with the result obtained by the previous method. We observe 
good agreement for large positive Lyapunov exponents, but increasing deviations for the 
smaller ones. These deviations can be eliminated by demanding higher numerical precision 
for the second method. 



IV. OTHER LATTICE FIELD THEORIES 

In this section, we want to apply our method to several other interacting fields. We 
calculate their largest Lyapunov exponents and investigate how they scale with energy, in 
order to determine the nature of the stochatic behavior of the various field theories. 
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A. U(l) Gauge Theory 



How characteristic is the observed behavior for the SU(2) gauge theory? Of course, one 
would expect a similar chaotic behavior for higher SU(N) gauge theories, because SU(2) is 
a subgroup of SU(N). However, one may expect a completely different behavior in the case 
of compact U(l), which has no self-interaction of gauge bosons in the continuum limit. The 
gauge group U(l) is easily obtained as a one-parameter abelian subgroup of SU(2). Not 
unexpectedly, we find quite different results in this case. For small values of the parameter 
g 2 Ea there is no discernable exponential divergence of adjacent field configurations. As 
shown in Figure 12, Aoa suddenly begins to grow rapidly for values g 2 Ea m 2. However, the 
dependence is highly nonlinear, reminiscent of the behavior of many nonlinear dynamical 
systems. This proves that the approximately linear relationship between Xqci and g 2 Ea 
found for the gauge group SU(2) is far from trivial. The latter is obviously associated with 
the nonabelian nature of this group. It is also worthwhile noting that Ao — > for a — > 
in U(l), i.e. the maximal Lyapunov exponent Ao for the compact U(l) lattice gauge theory 
vanishes in the continuum limit. 



B. SU(3) Gauge Theory 

Since SU(3) contains SU(2) as a subgroup, we will certainly expect that SU(3) gauge 
thoery is chaotic. This expectation is borne out by our study of the dynamics of SU(3) 
gauge fields on a lattice, which shows the same type of exponential instability 

0- While 

the basic procedure here is same as in SU(2), technically SU(3) is more complicated than 
SU(2). But nevertheless, since SU(3) is physically relevant as the gauge theory of the strong 
interaction, it is worthwhile to study its general dynamical properties. 

Whereas a SU(2) group element can be represented by a quaternion with only one redun- 
dant variable, there is no such simple representation for a general SU(3) group element. In 
principle we can use the exponential representation with 8 angles as independent variables. 
But as in the case of SU(2) this would involve many time consuming trigonometric manipu- 
lations required to carry out a group multiplication. In our simulations we have, therefore, 
directly used the matrix representation, i.e. we represent every group element by a unitary 
3x3 matrix. With the cost of 10 redundant variables for each link the group multipliaction is 
directly realized by matrix multiplication which only involves multiplications and additions. 

Another difference we shall point out is that it was necessary to be more careful in 
choosing the initial states here than in the case of SU(2), because the topological structure 
of the group space of SU(3) is more complicated. SU(3) is an eight-parameter Lie group 
and not all directions in group space are equivalent. Simply choosing angles at random in a 
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certain representation of the group parameters 13 lj can lead to a sampling of gauge fields that 
differs strongly from the thermal ensemble. It therefore proved necessary to choose the initial 



state by the heat-bath method 2j],|32j, which is straightforward but numerically expensive 
especially at low temperatures where the rejection rate is large. All other considerations, 
especially the scaling property (j2Tj) . carry over to SU(3). 

Our numerical results, obtained by evolving the thermalized initial condition, are pre- 
sented in Figure 13. They show again a close to linear dependence between Ao and g 2 E, but 
with a different slope than for SU(2). To good approximation we find: 

\ a « ±g 2 Ea. (40) 

We investigated several cases with non-thermalized initial conditions, specified by restricting 
some angles to a certain limited range. These initial conditions lead to exponential diver- 
gence rates lying above the line (|40|). We surmise that this reflects the lack of complete 
randomization of the field configuration during the limited time interval until the distance 
measure saturates. This suggests the importance of the thermalized initial condition, when 
the rescaling method is not applied. We have also tried a different definition of distance 

D E [U,U'} = i-£ I Etr^/^) -triU'J^) |, (41) 

iV P x i 

corresponding to the sum of the absolute value of local difference in the electric energy. The 
rise of In D E (t) is coincident with that exhibited by In D(t) except for the initial oscillatory 
region. This is within our expectation because the chaoticity of a system is an intrinsic 
property and it does not depend on a particular choice of the distance measure. 



C. Massless Scalar Field 

One question concerning the chaoticity of nonabelian gauge fields is whether the chaoticity 
is just a consequence of the nonlinearity of the field equations, or whether it is related to the 
particular form of the nonabelian gauge interaction. In this respect, we study the classical 
$ theory described by the Lagrangian 

L = dp&dPQ - /i 2 ^ - g 2 (&$) 2 . (42) 

The lattice formalism for interacting Higgs and gauge fields was given by Ambj0rn et al. 

We here study a special case, i.e. the massless limit of an iso-doublet complex scalar 
field, which has four real field components, without the gauge field. In the massless limit 
(/i = 0), just like in the case of a gauge theory, the corresponding classical lattice theory 
is scale invariant, in the sense that the self-coupling g 2 of the scalar field and the lattice 
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spacing a can be scaled out entirely and the system has no free parameter except the total 
energy. 

Using the same method as before, we measured Lyapunov exponents of the massless 
scalar field on the lattice. In Figure 14 we show X Q /g 2 E as function of g 2 Ea. Observe 
that the maximal Lyapunov exponent for the scalar field is much smaller than that for 
gauge fields with the same energy. Second, we find that the ratio Xo/g 2 E tends to zero at 
small g 2 Ea, suggesting a vanishing Lyapunov exponent in the continuum limit a — *■ 0. This 
result is consistent with our understanding of the relation between the maximal Lyapunov 
exponent and the damping rate of the long wavelength modes. The thermal perturbation 
theory calculation shows, unlike in the case of nonabelian gauge fields, that the damping 
rate of the long wavelength mode of the scalar field vanishes at the order of g 2 . This result 
shows that not all nonlinear classical continuum field theories are chaotic. 



D. Massive SU(2) Vector Fields 

This is the first step toward the understanding of chaoticity in the electroweak interaction, 
where the gauge symmetry is spontaneously broken by a Higgs field. We want to study how a 
mass term affects the chaotic behavior, in order to find a hint of how the interaction between 
the Higgs field and the gauge field may affect the chaoticity of the latter. In Figure 15, we 
show the largest Lyapunov exponents at various energies for theories with different vector 
boson mass. The same scaled variables are used as before, but here the mass parameter m 
cannot be scaled out. It is related to the mass M of the vector field in the continuum theory 
as M = mh/2a. For comparison, the results for the massless gauge theory are shown as solid 
squares, which are fit by a straight line with a coefficient of 4. The hollow squares are for 
m = 0.2, and the crosses are for m = 4. They more or less lie on a straight line. We observe 
that the effect of a mass term is to reduce the chaoticity of the system, or in other words, 
the mass term has a stabilizing effect on the trajectories of the field configurations. This 
effect is consistent with studies of the simple system (|12p. where the amount of chaoticity 



depends on the relative strength of the nonlinearity and the harmonic potential 11 1. 



E. Spontaneously Broken SU(2) Yang-Mills Theory 

The pure SU(2) and SU(3) Yang-Mills gauge theories are the primary tools for studying 
nonperturbative QCD related phenomena, such as properties of a hot quark-gluon plasma. 
The spontaneously broken SU(2) Yang-Mills theory, in which a charged scalar isodoublet 
field, the Higgs field, is coupled to the gauge boson, is a model of the electroweak gauge 
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theory. It is used to study the high-temperature limit of the standard model of particle 
physics, in particular, the phenomenon of baryon number violation at temperatures around 
the electroweak phase transition temperature of T = 200 GeV. 

In the limit of vanishing Weinberg angle, the Hamiltonian describing this model is given 

by 



H 



(fx 



\eiei + + $ f <i> + (A$)t(A$) - v 2 &$ + A($ t $) 2 , 

2 4 



where the dot symbol denotes time derivative, 




(43) 



(44) 



is a charged Higgs doublet in the fundamental representation of SU(2), and the gauge field is 
described by vector potentials A?(x) in the temporal axial gauge Aq = 0. In the numerical 
simulation we again discretize the gauge fields and represent them by link variables U X j, 
which are identified with real-valued quaternions. The Higgs doublet is also represented 
as a quaternion. However, its length (determinant) is also a dynamical degree of freedom, 
which we denote by R: 



R 2 = ^tr($ t $). 



(45) 



Following the notation of Amj0rn et al. |35l ] . we introduce a unit length quaternion V for the 
representation $ = RV. This factorization is a useful way to separate the gauge-invariant 
and the gauge-dependent degrees of freedom of the Higgs fields $. Nevertheless, in most 
formulas we shall use the matrices <£>, $t, for brevity. 

The lattice regularized version of the Hamiltonian (J43j) is based on the link- variables U x ^, 
defined as in eq. and the site- variables, $ x , or equivalently R X V X , for the rescaled Higgs 
field matrix 



(46) 



The Hamiltonian consists of three terms. The first one is the same as for the pure Yang-Mills 
theory 



.«w = — 



^ ]T tr (Ulfi x>i ) + E(l " ~tr U Xtij ) 

x,i x,ij 



(47) 



where U X) ij again stands for a plaquette quaternion. The second term in the Hamiltonian 
descibes the free Higgs doublet and its coupling to the gauge field 



H 



HW 



fa 

a 



(48) 
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where Higgs fields on neighboring sites x and x + i are coupled by the lattice transport 
operator (1 — U Xj i), that generates the gauge covariant derivative in the continuum limit. 
The last term in the Hamiltonian accounts for the self-interaction of the Higgs fields 

tf H = — £(^-l) 2 > (49) 

a X 

where R x is defined by eq. (|45|). The coupling parameters @q, /?h, and introduced into 
the lattice Hamiltonian are related to the parameters of the original continuum Hamiltonian 
as follows 

2 1 , Pr 2 4 

Ph a Pn Pg 

The vacuum expectation value of the Higgs field in the spontaneously broken phase is ob- 
tained from the minimum of the energy term Hn, and can be scaled to unity. The symmetry 
breaking, where (R x ) = 1, occurs when 

& > 3 " 27ft' (51) 



In this phase large Higgs and gauge boson masses are generated: 



Ml = 4h 2 \($ 2 ) = ^% (52) 



and 

The equations of motion derived from the lattice Hamiltonian 



M w = JftVW = (53) 



H = H W + H HW + H H , (54) 



are identical to ()24[ I25|) for the time derivatives of the components of the link variables, u x i 
and u xi . The three electric fields E c x i on each link are updated according to 

Ki = i £ tr [^(E^ - Dj >y )] + f^tr (W^^), ( 55 ) 
ag j z Pq z 

where the sum again runs over all four plaquettes attached to the link (x,i). Gauss' law 
now has the form 

DfE\ - ig(&T a <5> - &r a $) = 0, (56) 
and remains conserved under time evolution. Finally, the Higgs fields evolve according to 

= J2( U *,i®*+i + U x ,-^ x -i) - [6 + ^(Rl - l)]^, (57) 
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where the summation index i runs over all spatial directions (i = x,y,z). By introducing 
an auxiliary variable ^ x = $ x , we transform the set of second-order differential equations 
into a larger set of coupled first-order differential equations. 

Solving the evolution equations (|25I55I57|) we can study the dynamics of the coupled 
Yang-Mills-Higgs system, and see how the fields approach thermal equilibrium. In particular, 
we can investigate the divergence of trajectories in the configuration space. In addition to 
the already defined gauge-invariant distance measure (J2*9"|) for the distance between two 
Yang-Mills fields which is based on the magnetic energy, we measure the distance of Higgs 
field configurations by 

D B [$,&\ = ±Y,\ R *-K\> ( 58 ) 

X 

which is also gauge invariant. 

Depending on the relative strength of the gauge coupling g and Higgs self-coupling A 
different hierarchies in the thermalization rate can be obtained. Figure 16 shows the time 
evolution of the gauge field and Higgs-field distances (a) in a cross-coupled, A ~ g and (b) 
in a self-coupling dominated case, A g. The logarithmic growth rate, defined as 

h=j t lnD(t), (59) 

again depends on the average energy per plaquette. This dependence, counting only the 
energy per plaquette contained in the gauge field, is somewhat different in the above two 
cases, A ~ g and A g, as shown in Figure 17. At the upper end of the energy scale, 
h(E) is about the same for the two cases and also agrees with the value (J52*j) for the pure 
gauge field, h(E) is much smaller for the strongly coupled case (b) for low energies. We also 
note that the gauge field becomes chaotic faster than the Higgs field in case (b), while they 
become chaotic equally fast in case (a). 



V. THERMALIZATION OF GAUGE FIELDS 



A. Thermalization Time 

The universal exponential divergence of neighboring gauge field configurations for the 
gauge groups SU(2) and SU(3) implies that the entropy S of an ensemble of gauge fields 
grows linearly with time j^ ]: 

S(t) = S + t y £ Ai, (60) 

Ai>0 
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until the available microcanonical phase space is filled and the system is equilibrated. In 
perturbation theory, the average energy E per plaquette is related to the temperature as 

E « \{n 2 - 1)T [for SU(n)]. (61) 
It follows that the characteristic entropy growth rate, i.e. the thermalization rate is given 

by 

r,.^| Uls,TW " , 62 , 

[0Mg 2 T [SU(3)] 

where we have inserted our best-fit numerical values from eqs. (132)1 and (|4T)|) . Apart from 
a factor two, these values are remarkably close to those for the thermal damping rate for a 
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gauge boson at rest obtained by Braaten and Pisarski 

fifi o, N , J 0.175 g 2 T [SU(2)], 
7o = 6.635 -t—q T = < (63) 
247r y | 0.264 g 2 T [SU(3)]. 

At first sight, the relation r ~ 270 is quite surprising because these two quantities appear 
in totally different contexts and are calculated with different methods. On the one hand, the 
damping rate is the imaginary part of the self energy of a quasi-particle in a thermal gauge 
system and is calculated in the framework of effective quantum field theory. On the other 
hand, the Lyapunov exponent is a classical dynamical quantity describing the divergence of 
two classical gauge field trajectories. 

Although we do not yet know how to establish a direct relation between these two quan- 
tities, we understand that this similarity does not arise without reason. The two quantities, 
though very different from their contexts, both describe how fast a non-equilibrated gluon 
system approaches thermal equilibrium. The relevance of A n is clear from the above discus- 
sion. The connection of 7(0;) is apparent from the relation [3J| 

fM = ^- Y + c(co)e^\ (64) 

where f(ui,t) is time-dependent gauge boson distribution function, and the first term on 
the right-hand side is the Bose distribution. Obviously, 27 describes the rate of approach 
to equilibrium. It is possible to show that the gluon damping rate is basically a quantity of 
semi-classical origin 3c| . 

Finally, let us estimate the gluon thermalization quantitatively. Figure 18 shows r$ = Tq 1 
as function of temperature. For the gauge coupling constant g 2 we have used the renormal- 
ized running coupling constant of SU(3) gauge theory, which in the one-loop approximation 
is given by: 

» 2 < T > = Tu^JW' (65) 
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where A 200 MeV is the QCD scale parameter. For temperatures in the range T = 
300 — 500 MeV, which are realistically accessible in relativistic heavy-ion collisions, the 
thermalization time is less than 0.4 fm/c, or about 10~ 24 s. This promises the rapid formation 
of a locally thermalized gluon plasma in these collisions at sufficiently high energy. 

B. Self-Thermalization of Gauge Fields 

With the numerical tools in this context we can study the thermalization of gauge fields 
by direct simulation of their real-time evolution. Here we show that a nonabelian gauge 
field far off equilibrium approaches a thermally equilibrated state very rapidly. A thermally 
equilibrated state is the state in which the energy distribution over the microscopic degrees 
of freedom does not change with time and takes the thermal equilibrium form, which can 
be calculated from the canonical Gibbs ensemble. In order to study the process of thermal- 
ization, we must decide on which particular energy distribution we want to monitor. Here 
we choose as our monitor P(E M ), the distribution of magnetic energy on the elementary 
plaquettes over the whole lattice. The reason we prefer magnetic energy to electric energy is 
that the thermal distribution of the former can also be obtained using a heat bath algorithm, 
and so we can compare the results from time evolution and those from a heat bath method. 
P{Em) is obtained by counting the number of plaquettes that have magnetic energy be- 
tween E M — AE M /2 and E M + A.Em/2- ^E m should be sufficiently small so that P(E M ) 
is smooth but large enough to provide good counting statistics. Starting from an arbitrary 
intial state, we can measure the time evolution of P(E M ). If we find that it reaches some 
stable form we conclude that the system is thermally equilibrated. 

We numerically integrate the equations of motion for given initial conditions and measure 
the energy distribution functions P(E M ,t). The behavior in SU(2) and in SU(3) is rather 
similar. Here we concentrate on the SU(3) results. In Figure 19 we show the time evolution 
oiP(E M ) for SU(3). The corresponding initial averaged energy per plaquette is E = (E M ) = 
1.73. The plot shows \n(P(E M , t)/E M ) for reasons that will become clear below. The solid 
line denotes the initial distribution at t = 0. The dotted and short-dashed lines are for 
t = 0.5 and t — 1.5 respectively. The long-dashed line shows the final distribution reached 
at t — 3, whereafter no noticeable change is observed. The thermalization time is compatible 
with the inverse of the maximal Lyapunov exponent, which is Aq 1 = lO^E)^ 1 = 7.7 for 
SU(3). The final magnetic energy per plaquette is 0.84 which is almost half of the total 
energy. This is not always true on the lattice where the magnetic energy is limited from 
above because of the compact nature of the gauge group manifold. But in our example the 
total energy is small and the compactness is not a relevant influence. 

We have obtained equilibrated distributions P(E M ) for different initial energies in both 
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SU(2) and SU(3). These distributions can all be almost perfectly fit by 

P{E M ) = Uf(E M ) exp{-E M /T a ), (66) 

where M is a normalization constant, and T s is a "temperature" parameter. f(E) is the 
single plaquette phase space factor defined as 

f(E M )dE M = J S(E M - E[U])dE M dfi{U), (67) 

where E[U] = 2Re(n — trU) is the magnetic energy of a single free plaquette for SU(n) and 
dfi(U) is the Haar measure. Although the distribution looks like a Boltzmann distribution 
and T s like a temperature, we have shown that the real temperature of the system differs 
from T 8 . Their relation is shown in Figure 20 for SU(2). While the ratio T s /T — > 1 at high 
temperature, it decreases with T when T < 2 and finally it reaches | in the limit T = 0. 
We expect a similar behaviour for SU(3). This implies that in the low temperature limit 
the number of effective degrees of freedom is 2/3 of that in the high temperature limit. 



VI. CONCLUSION AND FUTURE DEVELOPMENTS 



The studies of lattice gauge fields reviewed here have considerably extended our knowledge 
of the properties of the dynamics of classical nonabelian gauge fields. This work began with 
the results of Matinyan, Savvidy and others more than a decade ago. We now know that the 
classical Yang- Mills fields are strongly chaotic at all energies. We have obtained the complete 
Lyapunov spectrum and shown that it scales even for very small lattice sizes. We are thus 
close to a full understanding of the classical Yang-Mills equations from the point of view 
of classical dynamics. We have also learned that, although other field theories show signs 
of chaotic behavior on the lattice, as well, their scaling properties imply that the chaoticity 
disappears in the continuum limit. This suggests that the chaotic behavior of nonabelian 
gauge theories is a nontrivial property of these theories. 

What are the consequences of the chaotic nature of the classical Yang-Mills theory? 
We have already studied its application to thermalization at finite temperature. Another 
possible application is the old idea 0| that the chaoticity may be related to the problem 
of color confinement. In order to investigate this conjecture, as well as other implications, 
more thoroughly, the following two possible directions of future research are important. 

In order to connect the results of classical gauge field dynamics to the underlying quantum 
field theory, we must understand the corrections from quantum effects. This means that 
we must develop appropriate semi-classical methods that link the classical limit to quantum 
physics. We have recently begun to explore a variational method using Gaussian wave 
packets for lattice gauge fields and found indications that the quantum corrections enhance 
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the chaoticity |38j. However, it is known that when a system is chaotic in the classical limit, 
the quantum evolution of a initially localized Gaussian wave packet tends to disperse the 
wave packet rapidly This is easily understood in the path integral formalulation of the 
quantum theory. A more realistic approach must account for this property. One possible 
approach is the Wigner function formulation mentioned in section 2.1. However, Wigner 
functions for gauge fields involve many subtleties jsjj, and it is not quite clear how the 
definition of a Wigner function should be best applied to gauge fields on a lattice. 

Secondly, the studies have so far been restricted to gauge fields and scalar fields. Since, 
in reality, the gauge fields are always coupled to fermion fields, a natural question is what 
role the fermions play in this context. Fermions cannot be treated classically. As a possible 
solution, on could consider a hybrid model decribing the interaction of the quantum evolution 
of matter fields with the classical evolution of the gauge fields. This requires the solution of 
the time-de pen dent Dirac equation on the lattice, which is equivalent to a quantum cellular 
automaton [40]. 

We must bear in mind that numerical studies of the real-time evolution of (quantum) field 
theories on lattices are in their infancy. The results obtained so far indicate that interesting 
results can be obtained from such investigations, and there is reason to hope that further 
progress is possible. 
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FIG. 1: Poincare surfaces of section of a single trajectory in the two-dimensional Yang-Mills 
model (z(t) = 0) for three values of the scaling parameter: r = 4.878 (left), r = 0.2 (center), and 
r = 0.0012 (right). 




FIG. 2: Distance distribution of points with the same distance D from a given SU(2) field config- 
uration for lattices of size 2 3 , 4 3 , and 20 3 . 
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FIG. 3: Evolution of the distance D(t) between neighboring random gauge field configurations 
for several average energies on a 20 3 lattice. The curves correspond, from top to bottom, to the 
parameters (a) 5 = 1,0.5,0.45,0.4,0.35; (b) S = 0.3,0.25,0.2,0.15,0.1. For every value of 5 two 
curves are shown, which are indistinguishable when 5 > 0.2. 
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FIG. 4: Initial oscillations in the distance between neighboring field configurations, before the 
exponential divergence sets in (for 5 = 0.03 and N = 10). 
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FIG. 5: Dependence of the exponential growth rate h on the average energy per plaquette E of the 
randomly chosen field configuration, for a = 0.5 and 4/g A = 1.1185, and for lattice sizes N = 20 
(open circles) and N = 6 (solid squares). The straight line through the origin is a least-squares fit. 
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FIG. 6: Scaling of the maximal Lyapunov exponent Ao with the variable g 2 Ea. At fixed value of 
the scaling variable, the exponent is independent of the value of g 2 over a wide range. The top 
part shows D(t) for different choices of g 2 , the bottom part shows the Lyapunov exponent. 




FIG. 7: Dependence of the evolution of the distance between field configurations on the lattice 
size N. The curves for N = 6, 10, 20, and 28 nearly coincide. All curves correspond to 5 = 0.4. 
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FIG. 8: The two largest Lyapunov exponents for SU(2) determined by the rescaling method. The 
average of the logarithmic scaling factors sf approaches the limit from above. 
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FIG. 9: Complete spectrum of 144 Lyapunov exponents for SU(2) gauge theory on a 2 3 lattice. The 
trajectories were followed up to time t/a = 200 (crosses) and t/a = 1000 (triangles). The central 
third fraction of Lyapunov exponents (enclosed between the vertical dashed lines) corresponds to 
the unphysical degrees of freedom that describe gauge transformations and deviations from Gauss' 
law. These exponents converge to zero in the limit t — > oo. 
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FIG. 10: Scaling of the Lyapunov spectrum with lattice size N. The solid line corresponds to a 3 3 
lattice; the dashed line is for a 2 3 lattice. Only the positive Lyapunov exponents are shown. The 
exponents Aj are scaled with the maximal Lyapunov exponent for each lattice size, and the index 
i is scaled with N 3 . 
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FIG. 11: Lyapunov spectrum for a 2 3 lattice obtained with the second scaling method working 
directly in the tangent space (hollow square), in comparison with the results based on the first 
method (solid triangle). 
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FIG. 12: Dependence of the scaled growth rate ha on the scaling parameter g 2 Ea for the gauge 
group U(l) on a 10 3 lattice. Note the highly nonlinear behavior and the rapid vanishing of ha in 
the limit g 2 Ea — > 0. 
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FIG. 13: Dependence of the scaled exponential growth rate ha = Aoa on the scaling parameter 
g 2 Ea for the gauge group SU(3). The calculations were performed on a 10 3 lattice. The dashed 
line depicts the fit by eq. (39). 



37 




FIG. 14: Maximal Lyapunov exponent for the massless <& 4 theory as function of the scaling variable 
g 2 Ea. 
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FIG. 15: Maximal Lyapunov exponent for the massive vector field theory as function of the scaling 
variable g 2 Ea for two different values of the vector boson mass: m = 0.2 (hollow squares), m = 4 
(crosses), and for comparison for m = (solid squares). 
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FIG. 16: Time evolution of the gauge field distance and the Higgs field distance for the sponta- 
neously broken SU(2) gauge theory. Part (a) shows a case where the gauge coupling g and the 
Higgs self-coupling A are about equal; part (b) corresponds to the case A > g. For the strongly 
coupled Higgs field (case b) most of the chaoticity resides in the gauge field. 
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FIG. 17: Maximal Lyapunov exponent of the Yang-Mills Higgs field as function of the variable 
g 2 Ea for the two cases discussed in the caption of Figure 16: (a) A = 0.9, (b) A = 9. The gauge 
coupling was g = 1.375, and the lattice spacing a = 0.5. The weakly coupled Higgs field follows 
very closely the results obtained for the pure gauge field. 
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FIG. 18: "Thermalization time" t° = A 1 for SU(3) as function of temperature T. The scale 
parameter was taken as A = 200 MeV. 
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FIG. 19: Time evolution of the distribution of magnetic plaquette energies Em for a SU(3) gauge 
field configuration. The plot shows the quantity P(Em)/E m on a logarithmic scale. The change in 
slope from the initial distribution (solid line) to the final distribution (long-dashed line) corresponds 
to the thermalization of all electric field modes, which has a cooling effect on the initially populated 
magnetic modes. 
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FIG. 20: Slope parameter ("apparent temperature") T s of the magnetic energy distribution P{Em) 
for the SU(2) gauge field as function of real temperature T. The plot shows the ratio T s /T, 
which rises from | to 1, because of contributions from the longitudinal plasma modes at higher 
temperature. 



42 



